Skip to main content
Version: Next

LightGBM - Quantile Regression for Drug Discovery (Scala)

Contents​

Overview & Background​

In pharmaceutical research and drug discovery, predicting the biological activity or potency of chemical compounds (Quantitative Structure-Activity Relationship, or QSAR) is a foundational task.

Traditional machine learning regression models optimize for Mean Squared Error (MSE), producing a single point estimate representing the conditional mean activity. However, in lead optimization and drug candidate selection, point estimates alone can be misleading:

  • Experimental assays have intrinsic measurement noise.
  • Novel chemical scaffolds often reside in sparse regions of chemical space (out-of-domain), where model confidence is naturally lower.
  • High variance and uncertainty can lead to costly laboratory synthesis and in vitro assay failures.

Quantile Regression addresses this challenge by estimating conditional percentiles (e.g., 20th percentile, 50th percentile / median, and 80th percentile) of the response distribution. Fitting models across multiple quantiles produces an uncertainty envelope (prediction interval) for every candidate compound. This empowers medicinal chemists to quantify risk, prioritize high-confidence candidates, and flag compounds requiring further experimental validation.

Runtime and dependency versions​

ComponentVersionNotes
Scala2.12.17Standalone build version; use the Scala 2.12 SynapseML artifact (synapseml_2.12)
Spark3.5.0Standalone build version; use a compatible Spark 3.5 runtime
SynapseML1.1.3Verify runtime supports this coordinate; managed runtimes may have different preinstalled versions
Hadoop connector (if using wasbs://)org.apache.hadoop:hadoop-azure:3.3.4Required only for standalone clusters reading wasbs:// blobs

Key Syntax Differences: PySpark vs. Spark Scala​

If you are transitioning from the Python SynapseML tutorial, keep these key differences in mind:

FeaturePython (PySpark)Scala (Spark)Explanation
Parameter ConfigurationLightGBMRegressor(alpha=0.5, objective="quantile")new LightGBMRegressor().setAlpha(0.5).setObjective("quantile")Scala uses the fluent setter pattern (.setParam()) instead of constructor keyword arguments.
Variable Immutabilitymodel = ...val model = ...Scala uses val for immutable bindings and var for mutable variables.
Array Definitions[0.8, 0.2]Array(0.8, 0.2)Scala uses typed collections (Array(...), Seq(...)).
Anonymous Functions[c for c in cols if c != "label"]cols.filter(_ != "label")Scala uses concise underscore _ syntax for lambdas.
Imports & Namespacesimport synapse.ml.lightgbm...import com.microsoft.azure.synapse.ml.lightgbm...Scala follows full JVM package hierarchy namespaces.

Step 1: Environment Setup and Dependencies​

To use LightGBM in Spark Scala, attach the SynapseML Maven coordinate to your Spark cluster or include it in your build configuration:

  • Maven Coordinate: com.microsoft.azure:synapseml_2.12:1.1.3
  • Spark Packages: com.microsoft.azure:synapseml_2.12:1.1.3
  • Repository: https://mmlspark.blob.core.windows.net/maven

Spark Shell / Databricks / Synapse Configuration​

When launching spark-shell or spark-submit, include the package:

spark-shell --packages com.microsoft.azure:synapseml_2.12:1.1.3 \
            --repositories https://mmlspark.blob.core.windows.net/maven

Note for Standalone Spark users (Option B — wasbs:// dataset): If you intend to use Option B (reading the public LibSVM dataset over Azure Blob Storage via wasbs://), you must also include the hadoop-azure connector. Managed cloud platforms (Databricks, Azure Synapse) pre-install this driver, but standalone Apache Spark does not include it by default. Add org.apache.hadoop:hadoop-azure:3.3.4 to --packages:

spark-shell \
  --packages com.microsoft.azure:synapseml_2.12:1.1.3,org.apache.hadoop:hadoop-azure:3.3.4 \
  --repositories https://mmlspark.blob.core.windows.net/maven

Without this, Spark will throw ClassNotFoundException: org.apache.hadoop.fs.azure.NativeAzureFileSystem$Secure.

Step 2: Spark Session and Imports​

In spark-shell or a managed Scala notebook, reuse the supplied spark session and import the classes below. The standalone application in Step 8 creates and stops its own session.

import org.apache.spark.sql.functions._
import org.apache.spark.ml.feature.VectorAssembler
import org.apache.spark.ml.evaluation.RegressionEvaluator
import com.microsoft.azure.synapse.ml.lightgbm.LightGBMRegressor

import spark.implicits._

Step 3: Dataset Preparation​

In QSAR modeling, compounds are typically represented by physicochemical descriptors or molecular fingerprints (e.g., Molecular Weight, LogP, Hydrogen Bond Donors/Acceptors, Topological Polar Surface Area, Rotatable Bonds) mapped to a biological potency target (such as pIC50=−log⁡10(IC50)pIC_{50} = -\log_{10}(IC_{50})).

Note: this tutorial uses a canonical target column name pIC50 across examples.

  • Option A (synthetic) uses pIC50.
  • Option B (LibSVM) renames the incoming label to pIC50 for consistency.
  • The standalone app uses pIC50 as the target column too.

Option A: Self-Contained Synthetic QSAR Dataset​

To run this tutorial immediately without external network dependencies, generate a synthetic QSAR dataset:

case class CompoundDescriptor(
  compound_id: String,
  molecular_weight: Double,
  logP: Double,
  hbd_count: Double,
  hba_count: Double,
  tpsa: Double,
  rotatable_bonds: Double,
  pIC50: Double // Bioactivity potency target
)

// Generate sample molecular descriptor data with heteroscedastic noise
val random = new scala.util.Random(42)
val sampleCompounds = (1 to 500).map { i =>
  val mw = 200.0 + random.nextDouble() * 350.0       // Molecular Weight (Da)
  val logP = -0.5 + random.nextDouble() * 5.5        // Octanol-water partition coefficient
  val hbd = random.nextInt(6).toDouble              // Hydrogen Bond Donors
  val hba = random.nextInt(10).toDouble             // Hydrogen Bond Acceptors
  val tpsa = 20.0 + random.nextDouble() * 120.0     // Topological Polar Surface Area (Ų)
  val rotBonds = random.nextInt(8).toDouble         // Rotatable Bonds

  // Synthetic QSAR response with scaffold-dependent variance (heteroscedasticity)
  val latentPotency = 4.0 + (0.005 * mw) + (0.4 * logP) - (0.15 * hbd) - (0.01 * tpsa)
  val noiseScale = 0.2 + 0.1 * (logP.abs)           // Uncertainty increases with extreme logP
  val noise = random.nextGaussian() * noiseScale
  val potency = latentPotency + noise

  CompoundDescriptor(s"CMPD-$i", mw, logP, hbd, hba, tpsa, rotBonds, potency)
}

val qsarDf = sampleCompounds.toDF()
qsarDf.show(5, truncate = false)

Option B: Public Triazines Benchmark Dataset (LibSVM)​

SynapseML also hosts the classic benchmark Triazines QSAR dataset (predicting inhibition of dihydrofolate reductase by pyrimidines).

LibSVM supplies a features vector and a label column (renamed to pIC50 below), so this path does not need VectorAssembler.

Standalone Apache Spark also needs Hadoop's Azure connector to read wasbs:// URLs. For Spark 3.5.0 with Hadoop 3.3.4, use --packages com.microsoft.azure:synapseml_2.12:1.1.3,org.apache.hadoop:hadoop-azure:3.3.4 with the Maven repository from Step 1. On managed clusters, use the connector supplied by the runtime or match the connector to the runtime's Hadoop version.

// Load benchmark Triazines QSAR dataset (requires cluster network connectivity)
val triazinesDf = spark.read
  .format("libsvm")
  .load("wasbs://publicwasb@mmlspark.blob.core.windows.net/triazines.scale.svmlight")

println(s"Total records in Triazines dataset: ${triazinesDf.count()}")
triazinesDf.printSchema()

// Rename LibSVM's default 'label' column to the canonical target column 'pIC50'
// so column names match the rest of this tutorial
val triazinesDfRenamed = triazinesDf.withColumnRenamed("label", "pIC50")

val Array(triazinesTrain, triazinesTest) = triazinesDfRenamed.randomSplit(Array(0.8, 0.2), seed = 1234L)
val triazinesModel = new LightGBMRegressor()
  .setObjective("quantile")
  .setAlpha(0.5)
  .setLabelCol("pIC50")
  .setFeaturesCol("features")
  .fit(triazinesTrain)

triazinesModel.transform(triazinesTest).select("pIC50", "prediction").show(5)

Tutorial Flow: The subsequent sections (Steps 4 through 7) follow Option A (qsarDf) to demonstrate how to perform custom feature engineering with VectorAssembler, multi-quantile uncertainty envelope modeling, and domain-specific bioactivity metric evaluation.

Step 4: Feature Assembly & Train/Test Split​

Assemble the molecular descriptor columns into a single Spark ML feature vector:

val featureCols = Array(
  "molecular_weight",
  "logP",
  "hbd_count",
  "hba_count",
  "tpsa",
  "rotatable_bonds"
)

val assembler = new VectorAssembler()
  .setInputCols(featureCols)
  .setOutputCol("features")

val assembledDf = assembler.transform(qsarDf)

// Split into training (80%) and testing (20%) datasets
val Array(trainData, testData) = assembledDf.randomSplit(Array(0.8, 0.2), seed = 1234L)
trainData.cache()
testData.cache()

println(s"Training set: ${trainData.count()} compounds")
println(s"Testing set:  ${testData.count()} compounds")

Step 5: Training Multi-Quantile LightGBM Models​

To construct an uncertainty envelope, train three separate LightGBMRegressor estimators configured with setObjective("quantile"):

  • α=0.20\alpha = 0.20: 20th percentile (conservative lower bound of compound potency).
  • α=0.50\alpha = 0.50: 50th percentile (median prediction, robust to outliers).
  • α=0.80\alpha = 0.80: 80th percentile (optimistic upper bound of compound potency).
// Helper method to create and configure a quantile regressor
def createQuantileRegressor(alpha: Double, predCol: String): LightGBMRegressor = {
  new LightGBMRegressor()
    .setObjective("quantile")
    .setAlpha(alpha)
    .setLabelCol("pIC50")
    .setFeaturesCol("features")
    .setPredictionCol(predCol)
    .setNumLeaves(31)
    .setNumIterations(100)
    .setLearningRate(0.05)
    .setMinDataInLeaf(10)
    .setSeed(42)
}

println("Training 20th percentile (Lower Bound) model...")
val modelQ20 = createQuantileRegressor(0.20, "pred_q20").fit(trainData)

println("Training 50th percentile (Median) model...")
val modelQ50 = createQuantileRegressor(0.50, "pred_q50_median").fit(trainData)

println("Training 80th percentile (Upper Bound) model...")
val modelQ80 = createQuantileRegressor(0.80, "pred_q80").fit(trainData)

Step 6: Generating the Uncertainty Envelope​

Transform the test data sequentially through all three models, then calculate the uncertainty width (q80−q20q_{80} - q_{20}):

val predictions = modelQ80.transform(
  modelQ50.transform(
    modelQ20.transform(testData)
  )
)

// Independent quantile fits do not guarantee monotonic ordering (q20 <= q50 <= q80).
// A crossing can reverse the endpoints or put the median outside them.
// See https://github.com/lightgbm-org/LightGBM/issues/3447.
val predictionsWithInterval = predictions
  .withColumn("is_crossed", $"pred_q20" > $"pred_q50_median" || $"pred_q50_median" > $"pred_q80")
  .withColumn("uncertainty_width", $"pred_q80" - $"pred_q20")
  .withColumn("within_interval", $"pIC50" >= $"pred_q20" && $"pIC50" <= $"pred_q80")
  .cache()

val numCrossedRows = predictionsWithInterval.filter($"is_crossed").count()
println(s"Rows with crossed quantiles: $numCrossedRows")

// Display sample predictions with uncertainty bounds and crossed flags
predictionsWithInterval
  .select("compound_id", "pIC50", "pred_q20", "pred_q50_median", "pred_q80", "uncertainty_width", "within_interval", "is_crossed")
  .show(10, truncate = false)

Interpretation for Medicinal Chemists​

For correctly ordered quantiles:

  • A small uncertainty_width means the estimated 20th and 80th percentiles are close. It does not establish model confidence or show that a compound is inside the training domain.
  • A large uncertainty_width means the estimated response interval is wide. These models do not separate assay noise from uncertainty in the fitted model.
  • A high pred_q20 is a high estimated lower response quantile, not a guaranteed minimum potency. Check held-out interval coverage and applicability to new compounds before using it for prioritization.

This synthetic example demonstrates the API, not a validated predictor of compound activity.

  • is_crossed flag: Compounds with reversed or crossed endpoints (is_crossed == true) have negative widths or inconsistent medians. They must not be treated as valid uncertainty intervals; taking absolute values or sorting does not establish advertised coverage.

Step 7: Model Evaluation & Validation​

Evaluate the median model using RMSE and MAE via RegressionEvaluator. Report coverage over all test rows and over the non-crossed subset separately. Neither is guaranteed to match the interval's nominal 60% coverage. The subset result describes only the retained rows, not the full test population.

// 1. Evaluate Median Model RMSE
val rmseEvaluator = new RegressionEvaluator()
  .setLabelCol("pIC50")
  .setPredictionCol("pred_q50_median")
  .setMetricName("rmse")

val rmse = rmseEvaluator.evaluate(predictionsWithInterval)
println(f"Median Model RMSE: $rmse%.4f")

// 2. Evaluate Median Model MAE
val maeEvaluator = new RegressionEvaluator()
  .setLabelCol("pIC50")
  .setPredictionCol("pred_q50_median")
  .setMetricName("mae")

val mae = maeEvaluator.evaluate(predictionsWithInterval)
println(f"Median Model MAE:  $mae%.4f")

// A positive width does not rule out a median outside the endpoints.
// Use the complete crossing diagnostic, not just uncertainty_width >= 0.
// See https://github.com/lightgbm-org/LightGBM/issues/3447.
val totalCount    = predictionsWithInterval.count()
require(totalCount > 0, "The test split must contain rows to evaluate coverage.")
val coverageCount = predictionsWithInterval.filter($"within_interval" === true).count()
val empiricalCoverage = (coverageCount.toDouble / totalCount.toDouble) * 100.0

val validRows          = predictionsWithInterval.filter(!$"is_crossed")
val validTotal         = validRows.count()
val validCoverageCount = validRows.filter($"within_interval" === true).count()
val validCoverage = if (validTotal > 0) {
  Some((validCoverageCount.toDouble / validTotal.toDouble) * 100.0)
} else {
  None
}

println(f"Rows with crossed quantiles        : $numCrossedRows (out of $totalCount)")
println(f"Empirical Coverage (all rows)      : $empiricalCoverage%.2f%%")
validCoverage match {
  case Some(coverage) =>
    println(f"Coverage on non-crossed rows      : $coverage%.2f%% ($validTotal rows)")
  case None =>
    println("Coverage on non-crossed rows      : not available, all rows have crossed quantiles")
}
println("The q20-q80 interval has nominal 60% coverage; measured coverage may differ.")

predictionsWithInterval.unpersist()
trainData.unpersist()
testData.unpersist()

Step 8: Standalone Spark Scala Application (spark-submit)​

To package this workflow into a standalone Scala application as requested in #731, create an sbt project. Put build.sbt in the project root and the application in src/main/scala/com/example/drugdiscovery/QSARQuantileApp.scala.

1. build.sbt​

name := "synapseml-lightgbm-qsar-standalone"
version := "1.0.0"
scalaVersion := "2.12.17"
resolvers += "SynapseML Maven Repo" at "https://mmlspark.blob.core.windows.net/maven"

val sparkVersion = "3.5.0"

libraryDependencies ++= Seq(
  "org.apache.spark" %% "spark-core" % sparkVersion % "provided",
  "org.apache.spark" %% "spark-sql"  % sparkVersion % "provided",
  "org.apache.spark" %% "spark-mllib" % sparkVersion % "provided",
  "com.microsoft.azure" % "synapseml_2.12" % "1.1.3"
)

2. Standalone Application (src/main/scala/com/example/drugdiscovery/QSARQuantileApp.scala)​

package com.example.drugdiscovery

import org.apache.spark.sql.SparkSession
import org.apache.spark.ml.feature.VectorAssembler
import org.apache.spark.ml.evaluation.RegressionEvaluator
import com.microsoft.azure.synapse.ml.lightgbm.LightGBMRegressor

object QSARQuantileApp {
  def main(args: Array[String]): Unit = {
    val spark = SparkSession.builder()
      .appName("QSAR-Quantile-Regression-Standalone")
      .getOrCreate()

    import spark.implicits._

    println("=== Running SynapseML LightGBM Quantile Regression Pipeline ===")

    // 1. Generate Synthetic Data
    val random = new scala.util.Random(42)
    val data = (1 to 1000).map { i =>
      val mw = 150.0 + random.nextDouble() * 400.0
      val logP = -1.0 + random.nextDouble() * 6.0
      val hbd = random.nextInt(5).toDouble
      val hba = random.nextInt(9).toDouble
      val pIC50 = 5.0 + 0.004 * mw + 0.35 * logP - 0.1 * hbd + random.nextGaussian() * 0.25
      (s"MOL_$i", mw, logP, hbd, hba, pIC50)
    }.toDF("id", "mw", "logP", "hbd", "hba", "pIC50")

    // 2. Assemble Features
    val assembler = new VectorAssembler()
      .setInputCols(Array("mw", "logP", "hbd", "hba"))
      .setOutputCol("features")

    val assembled = assembler.transform(data)
    val Array(train, test) = assembled.randomSplit(Array(0.8, 0.2), 42L)

    // 3. Train Quantile Models (10th, 50th, 90th percentiles for a nominal 80% prediction interval)
    val quantiles = Seq(
      (0.10, "pred_lower_10"),
      (0.50, "pred_median_50"),
      (0.90, "pred_upper_90")
    )

    var scoredTest = test
    for ((alpha, predCol) <- quantiles) {
      val model = new LightGBMRegressor()
        .setObjective("quantile")
        .setAlpha(alpha)
        .setLabelCol("pIC50")
        .setFeaturesCol("features")
        .setPredictionCol(predCol)
        .setNumLeaves(31)
        .setNumIterations(50)
        .setLearningRate(0.05)
        .fit(train)

      scoredTest = model.transform(scoredTest)
    }

    // 4. Compute Metrics
    val evaluator = new RegressionEvaluator()
      .setLabelCol("pIC50")
      .setPredictionCol("pred_median_50")
      .setMetricName("rmse")

    val rmse = evaluator.evaluate(scoredTest)
    println(f"Median Model RMSE: $rmse%.4f")

    scoredTest.select("id", "pIC50", "pred_lower_10", "pred_median_50", "pred_upper_90")
      .show(5, truncate = false)

    spark.stop()
  }
}

3. Execution via spark-submit​

Run these commands from the project root. Choose the submission command for your cluster manager. The standalone application uses synthetic data, so it does not need hadoop-azure. If you adapt it to read Option B's wasbs:// dataset, add the matching connector described in Step 1.

# Package the application
sbt package

# Run locally with an installed Spark 3.5 distribution
spark-submit \
  --class com.example.drugdiscovery.QSARQuantileApp \
  --master 'local[*]' \
  --deploy-mode client \
  --packages com.microsoft.azure:synapseml_2.12:1.1.3 \
  --repositories https://mmlspark.blob.core.windows.net/maven \
  target/scala-2.12/synapseml-lightgbm-qsar-standalone_2.12-1.0.0.jar

# Submit to a Spark standalone cluster; replace spark-master with your master host
spark-submit \
  --class com.example.drugdiscovery.QSARQuantileApp \
  --master spark://spark-master:7077 \
  --deploy-mode client \
  --packages com.microsoft.azure:synapseml_2.12:1.1.3 \
  --repositories https://mmlspark.blob.core.windows.net/maven \
  target/scala-2.12/synapseml-lightgbm-qsar-standalone_2.12-1.0.0.jar

# Submit to a configured YARN cluster with HADOOP_CONF_DIR or YARN_CONF_DIR set
spark-submit \
  --class com.example.drugdiscovery.QSARQuantileApp \
  --master yarn \
  --deploy-mode client \
  --packages com.microsoft.azure:synapseml_2.12:1.1.3 \
  --repositories https://mmlspark.blob.core.windows.net/maven \
  target/scala-2.12/synapseml-lightgbm-qsar-standalone_2.12-1.0.0.jar

For Databricks or Azure Synapse, use the platform's supported JAR submission mechanism and attach the SynapseML dependency. The commands above are local, Spark standalone, and YARN examples, not managed-platform submission commands.

Troubleshooting & common runtime errors​

  • ClassNotFoundException: org.apache.hadoop.fs.azure.NativeAzureFileSystem$Secure

    • Cause: missing Hadoop Azure connector on standalone Spark clusters when reading wasbs://.
    • Fix: add org.apache.hadoop:hadoop-azure:3.3.4 to --packages (or match your runtime's Hadoop version).
    • Example: 
      spark-shell \
        --packages com.microsoft.azure:synapseml_2.12:1.1.3,org.apache.hadoop:hadoop-azure:3.3.4 \
        --repositories https://mmlspark.blob.core.windows.net/maven

  • Quantile crossing (q20>q50q_{20} > q_{50} or q50>q80q_{50} > q_{80}):

    • Explanation: independent quantile fits can cross because each quantile regression model is trained separately without joint monotonic constraints. See maintainer explanation in LightGBM issue #3447.

Summary​

In this guide, you learned how to:

  1. Configure SynapseML LightGBM in Apache Spark Scala using current coordinates (1.1.3).
  2. Translate PySpark syntax to idiomatic Scala using fluent setter methods (.setParam()).
  3. Model biological activity (pIC50pIC_{50}) with Quantile Regression to estimate uncertainty intervals (q20,q50,q80q_{20}, q_{50}, q_{80}).
  4. Calculate empirical coverage and evaluate prediction accuracy with Spark ML's RegressionEvaluator.
  5. Package and submit a standalone Spark Scala LightGBM application using sbt and spark-submit.